
################################################################################
################################################################################
###                                                                          ###
### --- Filename: 3. Supplementary Material.R ------------------------------ ###
###                                                                          ###
### --- Description: This is the R script for the Supplementary Material for ###
### --- the paper "Economic Inequality and Public Support for -------------- ###
### --- Redistribution in Europe: A Cross-Sectional and Longitudinal ------- ###
### --- Multilevel Analysis" submitted to the IJPOR ------------------------ ###
###                                                                          ###
### --- Author: Valentin Velev --------------------------------------------- ###
###                                                                          ###
### --- Date: 11.12.2023 --------------------------------------------------- ###
###                                                                          ###
################################################################################
################################################################################

################################################################################
### --- 1. Preparation ----------------------------------------------------- ###
################################################################################

# install and load packages
required_packages <- c(
  # packages for loading and wrangling data
  "tidyverse", "haven", "readxl", "sjmisc", "sjlabelled", "fastDummies",
  # packages for (multiple) imputation
  "naniar", "imputeTS", "jomo", "mitml",
  # packages for statistical models
  "lme4", "lmerTest", "broom.mixed",
  # packages for post-estimation
  "influence.ME", "easystats", "sjPlot", "jtools", "interplot", "interactions",
  # other packages
  "ggpubr", "qqplotr", "NCmisc", "rstatix"
)

for (p in required_packages){
  
  if(!require(p, character.only=T)) install.packages(p, dependencies=T)
  library(p, character.only=T)
  
}

rm(p, required_packages)

# set seed to ensure replicability
set.seed(2946)

# global options
options(digits=6, scipen=999, max.print=10000)

# load necessary environment objects
load("Environment/data_with_NA.RData")
load("Environment/data.RData")
load("Environment/reg_dat.RData")
load("Environment/main_models.RData")
load("Environment/main_models_mi.RData")
load("Environment/additional_models_1.RData")
load("Environment/additional_models_2.RData")

################################################################################
### --- Table S1 ----------------------------------------------------------- ###
################################################################################

addmargins(table(data$country, data$year))

################################################################################
### --- Table S2 ----------------------------------------------------------- ###
################################################################################

as.data.frame(naniar::miss_var_summary(data_with_NA[-6:-9]))

################################################################################
### --- Table S3 ----------------------------------------------------------- ###
################################################################################

as.data.frame(rstatix::get_summary_stats(data[,c(-1:-3, -14:-20, -22:-26, -35:-37, -46:-49)], 
                                         show=c("n","mean","sd","min","max")))

################################################################################
### --- Table S4 ----------------------------------------------------------- ###
################################################################################

#
corr_mat <- 
  round(cor(data %>%
              subset(select=c("redist", "income", "lrscale", "employment", "edulvl",
                              "age", "male", "married", "badhealth", "unionmember",
                              "gini_disp_be", "gini_disp_we", "gdp_pc_be", "gdp_pc_we"))), 6)

#
corr_mat[upper.tri(corr_mat)] <- ""
corr_mat <- as.data.frame(corr_mat)
corr_mat

################################################################################
### --- Table S6 ----------------------------------------------------------- ###
################################################################################

# display results in regression table
tab_model(m3.1, m3.2, m3.3, m3.4, m3.5, m3.6,
          title="Model M3 - Additional country-level variables",
          dv.labels=c("M3.1", "M3.2", "M3.3", "M3.4", "M3.5", "M3.6"),
          string.est="b/SE", p.val="satterthwaite", p.style="stars", show.aic=T,
          show.se=T, show.ci=F, show.dev=T, collapse.se=T, digits=4, show.icc=F,
          digits.re=4, digits.rsq=4, file="Regression Tables/tab_s6.doc")

# display all variance and covariance components
data.frame(VarCorr(m3.1))
data.frame(VarCorr(m3.2))
data.frame(VarCorr(m3.3))
data.frame(VarCorr(m3.4))
data.frame(VarCorr(m3.5))
data.frame(VarCorr(m3.6))

# AIC & BIC values
jtools::summ(m3.1, digits=getOption("jtools-digits", default=2))
jtools::summ(m3.2, digits=getOption("jtools-digits", default=2))
jtools::summ(m3.3, digits=getOption("jtools-digits", default=2))
jtools::summ(m3.4, digits=getOption("jtools-digits", default=2))
jtools::summ(m3.5, digits=getOption("jtools-digits", default=2))
jtools::summ(m3.6, digits=getOption("jtools-digits", default=2))

################################################################################
### --- Table S7 ----------------------------------------------------------- ###
################################################################################

# display results in regression table
tab_model(m3_alt, m5_alt,
          title="Model M3 and M5 - Market Gini index", dv.labels=c("M3", "M5"),
          string.est="b/SE", p.val="satterthwaite", p.style="stars", show.aic=T,
          show.se=T, show.ci=F, show.dev=T, collapse.se=T, digits=4, show.icc=F,
          digits.re=4, digits.rsq=4, file="Regression Tables/tab_s7.doc")

# display all variance and covariance components
data.frame(VarCorr(m3_alt))
data.frame(VarCorr(m5_alt))

# AIC & BIC values
jtools::summ(m3_alt, digits=getOption("jtools-digits", default=2))
jtools::summ(m5_alt, digits=getOption("jtools-digits", default=2))

################################################################################
### --- Table S8 ----------------------------------------------------------- ###
################################################################################

save(reg_dat, m3, m5, file="models_for_outliertest.RData")

load("models_for_outliertest.RData")

# compute influence objects
influence_obj_m3 <- influence.ME::influence(m3, group="country") # note: takes around 30 min to compute
influence_obj_m5 <- influence.ME::influence(m5, group="country") # note: takes around 2 hrs to compute

# compute Cook's D
influence.ME::cooks.distance.estex(influence_obj_m3)
influence.ME::cooks.distance.estex(influence_obj_m5)

# compute DFBETAs
## Model M3
dfbetas_giniWE_m3 <- as.data.frame(dfbetas(influence_obj_m3)) %>%
  subset(select="gini_disp_we")

dfbetas_giniWE_m3

#dfbetas_giniWE_m3$country <- rownames(dfbetas_giniWE_m3)

## Model M5
dfbetas_giniWE_m5 <- as.data.frame(dfbetas(influence_obj_m5)) %>%
  subset(select="gini_disp_we")

dfbetas_giniWE_m5

# create Table S8
outlier_stats <-
  # create table with outlier statistics
  data.frame(Country = sort(unique(reg_dat$country)),
             Cooks_D_M3 = round(influence.ME::cooks.distance.estex(influence_obj_m3),4),
             Cooks_D_M5 = round(influence.ME::cooks.distance.estex(influence_obj_m5),4),
             DFBETAs_M3 = round(influence.ME::dfbetas.estex(influence_obj_m3,
                                                            parameters="gini_disp_we"),4),
             DFBETAs_M5 = round(influence.ME::dfbetas.estex(influence_obj_m5,
                                                            parameters="gini_disp_we"),4)) %>%
  # add stars (cut-off values)
  mutate(Cooks_D_cutoff_M3 = case_when(Cooks_D_M3 >= 0.1379 ~ "*",
                                       T ~ " "), .after=Cooks_D_M3) %>%
  mutate(Cooks_D_cutoff_M5 = case_when(Cooks_D_M5 >= 0.1379 ~ "*",
                                       T ~ " "), .after=Cooks_D_M5) %>%
  mutate(DFBETAs_cutoff_M3 = case_when(DFBETAs_M3 <= -0.3714 | 
                                         DFBETAs_M3 >= 0.3714 ~ "*",
                                       T ~ " "), .after=DFBETAs_M3) %>%
  mutate(DFBETAs_cutoff_M5 = case_when(DFBETAs_M5 <= -0.3714 | 
                                         DFBETAs_M5 >= 0.3714 ~ "*",
                                       T ~ " "), .after=DFBETAs_M5) %>%
  # paste statistic and stars together
  unite("cooksd_m3", c(Cooks_D_M3, Cooks_D_cutoff_M3), sep="") %>%
  unite("cooksd_m5", c(Cooks_D_M5, Cooks_D_cutoff_M5), sep="") %>%
  unite("dfbetas_m3", c(DFBETAs_M3, DFBETAs_cutoff_M3), sep="") %>%
  unite("dfbetas_m5", c(DFBETAs_M5, DFBETAs_cutoff_M5), sep="")

# display table
outlier_stats

################################################################################
### --- Table S9 ----------------------------------------------------------- ###
################################################################################

# re-estimating Model M3 (Table 1 in paper) without outliers
## excluding BG
m3_no_BG <- lmer(redist ~
                   income_cgm+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
                   married+badhealth+unionmember+
                   gini_disp_be+gini_disp_we+gdp_pc_be+gdp_pc_we+
                   year+(1|country)+(1|country:year),
                 data=reg_dat[reg_dat$country!="BG",], REML=T,
                 control=lmerControl(optimizer="bobyqa",
                                     optCtrl=list(maxfun=1e6),
                                     calc.derivs=F))

## excluding DE
m3_no_DE <- lmer(redist ~
                   income_cgm+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
                   married+badhealth+unionmember+
                   gini_disp_be+gini_disp_we+gdp_pc_be+gdp_pc_we+
                   year+(1|country)+(1|country:year),
                 data=reg_dat[reg_dat$country!="DE",], REML=T,
                 control=lmerControl(optimizer="bobyqa",
                                     optCtrl=list(maxfun=1e6),
                                     calc.derivs=F))

## excluding NO
m3_no_NO <- lmer(redist ~
                   income_cgm+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
                   married+badhealth+unionmember+
                   gini_disp_be+gini_disp_we+gdp_pc_be+gdp_pc_we+
                   year+(1|country)+(1|country:year),
                 data=reg_dat[reg_dat$country!="NO",], REML=T,
                 control=lmerControl(optimizer="bobyqa",
                                     optCtrl=list(maxfun=1e6),
                                     calc.derivs=F))

## excluding UK
m3_no_UK <- lmer(redist ~
                   income_cgm+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
                   married+badhealth+unionmember+
                   gini_disp_be+gini_disp_we+gdp_pc_be+gdp_pc_we+
                   year+(1|country)+(1|country:year),
                 data=reg_dat[reg_dat$country!="UK",], REML=T,
                 control=lmerControl(optimizer="bobyqa",
                                     optCtrl=list(maxfun=1e6),
                                     calc.derivs=F))

## excluding all outliers (BG, DE, NO, UK)
m3_no_outliers <- lmer(redist ~
                         income_cgm+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
                         married+badhealth+unionmember+
                         gini_disp_be+gini_disp_we+gdp_pc_be+gdp_pc_we+
                         year+(1|country)+(1|country:year),
                       data=reg_dat[reg_dat$country!="BG" &
                                      reg_dat$country!="DE" &
                                      reg_dat$country!="NO" &
                                      reg_dat$country!="UK",], REML=T,
                       control=lmerControl(optimizer="bobyqa",
                                           optCtrl=list(maxfun=1e6),
                                           calc.derivs=F))

# display results in regression table (Table S9 in Supplementary Material)
tab_model(m3_no_BG, m3_no_DE, m3_no_NO, m3_no_UK, m3_no_outliers,
          title="Model M3 - No outliers",
          dv.labels=c("M3 w/o BG", "M3 w/o DE", "M3 w/o NO", "M3 w/o UK", "M3 w/o outliers"),
          string.est="b/SE", p.val="satterthwaite", p.style="stars", show.aic=T,
          show.se=T, show.ci=F, show.dev=T, collapse.se=T, digits=4, show.icc=F,
          digits.re=4, digits.rsq=4, file="Regression Tables/tab_s9.doc")

# display all variance and covariance components
data.frame(VarCorr(m3_no_BG))
data.frame(VarCorr(m3_no_DE))
data.frame(VarCorr(m3_no_NO))
data.frame(VarCorr(m3_no_UK))
data.frame(VarCorr(m3_no_outliers))

# AIC & BIC values
jtools::summ(m3_no_BG, digits=getOption("jtools-digits", default=2))
jtools::summ(m3_no_DE, digits=getOption("jtools-digits", default=2))
jtools::summ(m3_no_NO, digits=getOption("jtools-digits", default=2))
jtools::summ(m3_no_UK, digits=getOption("jtools-digits", default=2))
jtools::summ(m3_no_outliers, digits=getOption("jtools-digits", default=2))

################################################################################
### --- Table S10 ----------------------------------------------------------- ###
################################################################################

# re-estimating Model M5 (Table 1 in paper) without outliers
## excluding BG
m5_no_BG <- lmer(redist ~
                   income_cwc+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
                   married+badhealth+unionmember+
                   gini_disp_be+gini_disp_we+gdp_pc_be+gdp_pc_we+
                   (gini_disp_be*income_cwc)+(gini_disp_we*income_cwc)+
                   year+(1+income_cwc|country)+(1+income_cwc|country:year),
                 data=reg_dat[reg_dat$country!="BG",], REML=T,
                 control=lmerControl(optimizer="bobyqa",
                                     optCtrl=list(maxfun=1e6),
                                     calc.derivs=F))

## excluding DE
m5_no_DE <- lmer(redist ~
                   income_cwc+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
                   married+badhealth+unionmember+
                   gini_disp_be+gini_disp_we+gdp_pc_be+gdp_pc_we+
                   (gini_disp_be*income_cwc)+(gini_disp_we*income_cwc)+
                   year+(1+income_cwc|country)+(1+income_cwc|country:year),
                 data=reg_dat[reg_dat$country!="DE",], REML=T,
                 control=lmerControl(optimizer="bobyqa",
                                     optCtrl=list(maxfun=1e6),
                                     calc.derivs=F))

## excluding NO
m5_no_NO <- lmer(redist ~
                   income_cwc+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
                   married+badhealth+unionmember+
                   gini_disp_be+gini_disp_we+gdp_pc_be+gdp_pc_we+
                   (gini_disp_be*income_cwc)+(gini_disp_we*income_cwc)+
                   year+(1+income_cwc|country)+(1+income_cwc|country:year),
                 data=reg_dat[reg_dat$country!="NO",], REML=T,
                 control=lmerControl(optimizer="bobyqa",
                                     optCtrl=list(maxfun=1e6),
                                     calc.derivs=F))

## excluding UK
m5_no_UK <- lmer(redist ~
                   income_cwc+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
                   married+badhealth+unionmember+
                   gini_disp_be+gini_disp_we+gdp_pc_be+gdp_pc_we+
                   (gini_disp_be*income_cwc)+(gini_disp_we*income_cwc)+
                   year+(1+income_cwc|country)+(1+income_cwc|country:year),
                 data=reg_dat[reg_dat$country!="UK",], REML=T,
                 control=lmerControl(optimizer="bobyqa",
                                     optCtrl=list(maxfun=1e6),
                                     calc.derivs=F))

## excluding all outliers (BG, DE, NO, UK)
m5_no_outliers <- lmer(redist ~
                         income_cwc+lrscale_cgm+age_cgm+I(age_cgm^2)+male+edulvl+employment+
                         married+badhealth+unionmember+
                         gini_disp_be+gini_disp_we+gdp_pc_be+gdp_pc_we+
                         (gini_disp_be*income_cwc)+(gini_disp_we*income_cwc)+
                         year+(1+income_cwc|country)+(1+income_cwc|country:year),
                       data=reg_dat[reg_dat$country!="BG" &
                                      reg_dat$country!="DE" &
                                      reg_dat$country!="NO" &
                                      reg_dat$country!="UK",], REML=T,
                       control=lmerControl(optimizer="bobyqa",
                                           optCtrl=list(maxfun=1e6),
                                           calc.derivs=F))

# display results in regression table (Table S10 in Supplementary Material)
tab_model(m5_no_BG, m5_no_DE, m5_no_NO, m5_no_UK, m5_no_outliers,
          title="Model M5 - No outliers",
          dv.labels=c("M5 w/o BG", "M5 w/o DE", "M5 w/o NO", "M5 w/o UK", "M5 w/o outliers"),
          string.est="b/SE", p.val="satterthwaite", p.style="stars", show.aic=T,
          show.se=T, show.ci=F, show.dev=T, collapse.se=T, digits=4, show.icc=F,
          digits.re=4, digits.rsq=4, file="Regression Tables/tab_s10.doc")

# display all variance and covariance components
data.frame(VarCorr(m5_no_BG))
data.frame(VarCorr(m5_no_DE))
data.frame(VarCorr(m5_no_NO))
data.frame(VarCorr(m5_no_UK))
data.frame(VarCorr(m5_no_outliers))

# AIC & BIC values
jtools::summ(m5_no_BG, digits=getOption("jtools-digits", default=2))
jtools::summ(m5_no_DE, digits=getOption("jtools-digits", default=2))
jtools::summ(m5_no_NO, digits=getOption("jtools-digits", default=2))
jtools::summ(m5_no_UK, digits=getOption("jtools-digits", default=2))
jtools::summ(m5_no_outliers, digits=getOption("jtools-digits", default=2))

################################################################################
### --- Supplementary Material - Figures ----------------------------------- ###
################################################################################

################################################################################
### --- Figure S1 ---------------------------------------------------------- ###
################################################################################

#
coefs_m3 <- 
  broom.mixed::tidy(m3, conf.int=T, conf.level=.95)[,c("term", "estimate", "conf.low", "conf.high")] %>%
  sjmisc::rename_variables(term=Variable, estimate=Estimate, conf.low=ci95_low_m3, conf.high=ci95_high_m3)

coefs_m3 <- coefs_m3[2:27,]

#
coefs_m3_mi <- m3_mi %>% 
  rownames_to_column("Variable") %>%
  mutate(ci95_low_m3 = Estimate - 1.96 * Std.Error, ci95_high_m3 = Estimate + 1.96 * Std.Error)

coefs_m3_mi <- coefs_m3_mi[2:27,c(1,2,9,10)]

#
coef_plot1_df <- rbind(coefs_m3, coefs_m3_mi) %>% 
  mutate(model = ifelse(row_number() <= 26, "LWD", "MI"))

# create Figure S1
figure_s1 <- ggplot(data=coef_plot1_df[-c(19:26,45:52),],
                    aes(x=Estimate, y=Variable, xmin=ci95_low_m3, xmax=ci95_high_m3, color=model)) +
  geom_vline(xintercept=0, color="black", linetype="dashed", linewidth=.3) +
  geom_pointrange(position=position_dodge(width=-.9), size=.1) +
  scale_y_discrete(limits=c("gdp_pc_we", "gdp_pc_be", "gini_disp_we", "gini_disp_be", 
                            "unionmember1", "badhealth1", "married1", "male1", "I(age_cgm^2)", 
                            "age_cgm", "edulvl5", "edulvl4", "edulvl3", "edulvl2", "employment2", 
                            "employment1", "lrscale_cgm", "income_cgm"),
                   labels=c("unionmember1"="Union member",
                            "married1"="Married",
                            "male1"="Male",
                            "lrscale_cgm"="L-R self-placement",
                            "income_cgm"="Household income",
                            "I(age_cgm^2)"="Age squared",
                            "gini_disp_we"="Gini Index [WE]",
                            "gini_disp_be"="Gini Index [BE]",
                            "gdp_pc_we"="GDP/capita [WE]",
                            "gdp_pc_be"="GDP/capita [BE]",
                            "employment2"="Employment: unemployed",
                            "employment1"="Employment: employed",
                            "edulvl5"="Education: ISCED 5-6",
                            "edulvl4"="Education: ISCED 4",
                            "edulvl3"="Education: ISCED 3",
                            "edulvl2"="Education: ISCED 2",
                            "badhealth1"="Bad health",
                            "age_cgm"="Age")) +
  scale_color_manual(values=c("LWD"="black", "MI"="darkgrey"), name="Missing Data \nTreatment") +
  labs(x="Coefficient (\u03b2)", y="Variable") +
  theme(axis.title.x = element_text(size=9),
        axis.title.y = element_text(size=9),
        axis.text = element_text(size=6),
        plot.background = element_blank(),
        panel.background = element_rect(fill="white", colour="black"),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        legend.title = element_text(size=9),
        legend.text = element_text(size=7),
        legend.key = element_rect(fill="white"),
        text = element_text(family="Serif", face="bold"))

# save plot as .png file
ggsave(filename="Figures/figure_s1.png", plot=figure_s1, width=6, height=4, units="in", dpi=1500)

################################################################################
### --- Figure S2 ---------------------------------------------------------- ###
################################################################################

#
coefs_m5 <- 
  broom.mixed::tidy(m5, conf.int=T, conf.level=.95)[,c("term", "estimate", "conf.low", "conf.high")] %>%
  sjmisc::rename_variables(term=Variable, estimate=Estimate, conf.low=ci95_low_m5, conf.high=ci95_high_m5)

coefs_m5 <- coefs_m5[2:29,]

#
coefs_m5_mi <- m5_mi %>%
  rownames_to_column("Variable") %>%
  mutate(ci95_low_m5 = Estimate - 1.96 * Std.Error, ci95_high_m5 = Estimate + 1.96 * Std.Error)

coefs_m5_mi <- coefs_m5_mi[2:29,c(1,2,9,10)]

#
coef_plot2_df <- rbind(coefs_m5, coefs_m5_mi) %>% 
  mutate(model = ifelse(row_number() <= 28, "LWD", "MI"))

# create Figure S2
figure_s2 <- ggplot(data=coef_plot2_df[-c(19:26,47:54),], 
                    aes(x=Estimate, y=Variable, xmin=ci95_low_m5, xmax=ci95_high_m5, color=model)) +
  geom_vline(xintercept=0, color="black", linetype="dashed", linewidth=.3) +
  geom_pointrange(position=position_dodge(width=-.9), size=.1) +
  scale_y_discrete(limits=c("income_cwc:gini_disp_we", "income_cwc:gini_disp_be", "gdp_pc_we", 
                            "gdp_pc_be", "gini_disp_we", "gini_disp_be", "unionmember1", 
                            "badhealth1", "married1", "male1", "I(age_cgm^2)", "age_cgm", 
                            "edulvl5", "edulvl4", "edulvl3", "edulvl2", "employment2", 
                            "employment1", "lrscale_cgm", "income_cwc"),
                   labels=c("unionmember1"="Union member",
                            "married1"="Married",
                            "male1"="Male",
                            "lrscale_cgm"="L-R self-placement",
                            "income_cwc"="Household income",
                            "I(age_cgm^2)"="Age squared",
                            "gini_disp_we"="Gini Index [WE]",
                            "gini_disp_be"="Gini Index [BE]",
                            "gdp_pc_we"="GDP/capita [WE]",
                            "gdp_pc_be"="GDP/capita [BE]",
                            "income_cwc:gini_disp_we"="Income x Gini Index [WE]", 
                            "income_cwc:gini_disp_be"="Income x Gini Index [BE]",
                            "employment2"="Employment: unemployed",
                            "employment1"="Employment: employed",
                            "edulvl5"="Education: ISCED 5-6",
                            "edulvl4"="Education: ISCED 4",
                            "edulvl3"="Education: ISCED 3",
                            "edulvl2"="Education: ISCED 2",
                            "badhealth1"="Bad health",
                            "age_cgm"="Age")) +
  scale_color_manual(values=c("LWD"="black", "MI"="darkgrey"), name="Missing Data \nTreatment") +
  labs(x="Coefficient (\u03b2)", y="Variable") +
  theme(axis.title.x = element_text(size=9),
        axis.title.y = element_text(size=9),
        axis.text = element_text(size=6),
        plot.background = element_blank(),
        panel.background = element_rect(fill="white", colour="black"),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        legend.title = element_text(size=9),
        legend.text = element_text(size=7),
        legend.key = element_rect(fill="white"),
        text = element_text(family="Serif", face="bold"))

# save plot as .png file
ggsave(filename="Figures/figure_s2.png", plot=figure_s2, width=6, height=4, units="in", dpi=1500)

################################################################################
### --- Figure S3 ---------------------------------------------------------- ###
################################################################################

# create Figure S3
fig_s3 <- ggplot(data=data, aes(x=gini_disp_be, y=redist_be)) +
  geom_smooth(data=data[data$country=="AT",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="BE",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="BG",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="CH",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="CY",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="CZ",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="DE",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="DK",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="EE",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="ES",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="FI",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="FR",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="GR",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="HR",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="HU",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="IE",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="IS",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="IT",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="LT",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="LV",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="NL",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="NO",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="PL",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="PT",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="RO",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="SE",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="SI",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="SK",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(data=data[data$country=="UK",], aes(x=gini_disp, y=redist_C), 
              method="lm", se=F, color="grey", alpha=.35, linewidth=.9, inherit.aes=F) +
  geom_smooth(formula=y~x, method="lm", color="black", se=F) +
  geom_point(size=1.75, color="black") +
  geom_text(data=data, aes(label=country), hjust=1.3, vjust=-.35, size=2.45, family="Serif") +
  labs(x="Income Inequality", 
       y="Demand for Redistribution") +
  theme(plot.background = element_blank(),
        panel.background = element_rect(fill="white", colour="black"),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        text = element_text(family="Serif", face="bold"))

# save plot as .png file
ggsave(filename="Figures/fig_s3.png", plot=fig_s3, width=6, height=4, units="in", dpi=1500)

################################################################################
### --- Figure S4 ---------------------------------------------------------- ###
################################################################################

# aggregate data
data_agg <- data %>%
  mutate(year = as.integer(year)) %>%
  group_by(country, year) %>%
  summarise(across(c(redist, gini_disp), \(x) mean(x, na.rm=T))) %>%
  ungroup()

# create Figure S4
fig_s4 <- ggplot(data=data_agg, aes(x=year)) +
  geom_line(aes(y=redist/(max(redist)/max(gini_disp)), linetype="line1"), linewidth=.25) +
  geom_point(aes(y=redist/(max(redist)/max(gini_disp))), size=.75, show.legend=F, color="#333333") +
  geom_line(aes(y=gini_disp, linetype="line2"), linewidth=.25) +
  geom_point(aes(y=gini_disp), size=.75, show.legend=F, color="darkgrey") +
  scale_y_continuous(name="Gini Index", 
                     sec.axis = sec_axis(trans=~.*(max(data_agg$redist)/max(data_agg$gini_disp)), 
                                         name="Demand for Redistribution")) +
  scale_linetype_manual(name="", 
                        values=c("line1"=1,"line2"=2), 
                        labels=c("Demand for Redistribution","Gini Index")) +
  facet_wrap(~country, scales="free_y") +
  labs(x="Year") +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12),
        axis.text = element_text(size=7),
        strip.text.x = element_text(size=7),
        strip.background = element_blank(),
        plot.background = element_blank(),
        panel.background = element_rect(fill="white", colour="black"),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        text = element_text(family="Serif", face="bold"),
        legend.position = "bottom",
        legend.background = element_rect(fill="white", colour="black"),
        legend.key = element_blank())

# save plot as .png file
ggsave(filename="Figures/fig_s4.png", plot=fig_s4, width=10, height=8, units="in", dpi=1500)

################################################################################
### --- Figure S5 ---------------------------------------------------------- ###
################################################################################

#
jn_plot1 <- johnson_neyman(m5, pred="income_cwc", modx="gini_disp_be", sig.color="black",
                           insig.color="grey", line.thickness=.75, title=NULL, control.fdr=T)[["plot"]] +
  labs(x="Gini Index [BE]", y="CME of Household Income") +
  theme(axis.title.x = element_text(size=9),
        axis.title.y = element_text(size=9),
        axis.text = element_text(size=7),
        legend.text = element_text(size=8),
        text = element_text(family="Serif", face="bold"))

#
jn_plot2 <- johnson_neyman(m5, pred="income_cwc", modx="gini_disp_we", sig.color="black",
                           insig.color="grey", line.thickness=.75, title=NULL, control.fdr=T)[["plot"]] +
  labs(x="Gini Index WE]", y="CME of Household Income") +
  theme(axis.title.x = element_text(size=9),
        axis.title.y = element_text(size=9),
        axis.text = element_text(size=7),
        legend.text = element_text(size=8),
        text = element_text(family="Serif", face="bold"))

# append plots (create Figure S5)
fig_s5 <- ggarrange(jn_plot1, jn_plot2, ncol=2, common.legend=T, legend="bottom")

# save plot .png file
ggsave(filename="Figures/fig_s5.png", plot=fig_s5, width=8, height=3.5, units="in", dpi=1750)

################################################################################
### --- Figure S6 ---------------------------------------------------------- ###
################################################################################

# extract upper-level residuals from Model M3
resid_m3_country <- ranef(m3)[["country"]]

resid_m3_country_year <- ranef(m3)[["country:year"]]

# create country-level P-P plots
pp_plot_m3_country <- ggplot(data=resid_m3_country,
                             aes(sample=resid_m3_country$`(Intercept)`)) +
  qqplotr::stat_pp_point(colour="#333333", shape=16) +
  qqplotr::stat_pp_line(colour="black") +
  labs(x="Theoretical Cumulative Distribution",
       y="Empirical Cumulative Distribution") +
  theme(axis.title.x = element_text(size=8),
        axis.title.y = element_text(size=8),
        axis.text = element_text(size=6),
        plot.background = element_blank(),
        panel.background = element_rect(fill="white", colour="black"),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        text = element_text(family="Serif", face="bold"))

# create country-year-level P-P plots
pp_plot_m3_country_year <- ggplot(data=resid_m3_country_year,
                                  aes(sample=resid_m3_country_year$`(Intercept)`)) +
  qqplotr::stat_pp_point(colour="#333333", shape=16) +
  qqplotr::stat_pp_line(colour="black") +
  labs(x="Theoretical Cumulative Distribution",
       y="Empirical Cumulative Distribution") +
  theme(axis.title.x = element_text(size=8),
        axis.title.y = element_text(size=8),
        axis.text = element_text(size=6),
        plot.background = element_blank(),
        panel.background = element_rect(fill="white", colour="black"),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        text = element_text(family="Serif", face="bold"))

# append plots (create Figure S6)
fig_s6 <- ggarrange(pp_plot_m3_country, pp_plot_m3_country_year,
                    ncol=2, labels=c("a","b"), 
                    font.label=list(size=12, color="black"))

# save plot .png file
ggsave(filename="Figures/fig_s6.png", plot=fig_s6, width=6, height=4, units="in", dpi=1500)

################################################################################
### --- Figure S7 ---------------------------------------------------------- ###
################################################################################

# extract upper-level residuals from Model M5
resid_m5_country <- ranef(m5)[["country"]]

resid_m5_country_year <- ranef(m5)[["country:year"]]

# create country-level P-P plots
pp_plot_m5_country1 <- ggplot(data=resid_m5_country,
                              aes(sample=resid_m5_country$`(Intercept)`)) +
  qqplotr::stat_pp_point(colour="#333333", shape=16) +
  qqplotr::stat_pp_line(colour="black") +
  labs(x="Theoretical Cumulative Distribution",
       y="Empirical Cumulative Distribution") +
  theme(axis.title.x = element_text(size=6),
        axis.title.y = element_text(size=6),
        axis.text = element_text(size=5),
        plot.background = element_blank(),
        panel.background = element_rect(fill="white", colour="black"),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        text = element_text(family="Serif", face="bold"))

pp_plot_m5_country2 <- ggplot(data=resid_m5_country,
                              aes(sample=resid_m5_country$income_cwc)) +
  qqplotr::stat_pp_point(colour="#333333", shape=16) +
  qqplotr::stat_pp_line(colour="black") +
  labs(x="Theoretical Cumulative Distribution",
       y="Empirical Cumulative Distribution") +
  theme(axis.title.x = element_text(size=6),
        axis.title.y = element_text(size=6),
        axis.text = element_text(size=5),
        plot.background = element_blank(),
        panel.background = element_rect(fill="white", colour="black"),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        text = element_text(family="Serif", face="bold"))

# create country-year-level P-P plots
pp_plot_m5_country_year1 <- ggplot(data=resid_m5_country_year,
                                   aes(sample=resid_m5_country_year$`(Intercept)`)) +
  qqplotr::stat_pp_point(colour="#333333", shape=16) +
  qqplotr::stat_pp_line(colour="black") +
  labs(x="Theoretical Cumulative Distribution",
       y="Empirical Cumulative Distribution") +
  theme(axis.title.x = element_text(size=6),
        axis.title.y = element_text(size=6),
        axis.text = element_text(size=5),
        plot.background = element_blank(),
        panel.background = element_rect(fill="white", colour="black"),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        text = element_text(family="Serif", face="bold"))

pp_plot_m5_country_year2 <- ggplot(data=resid_m5_country_year,
                                   aes(sample=resid_m5_country_year$income_cwc)) +
  qqplotr::stat_pp_point(colour="#333333", shape=16) +
  qqplotr::stat_pp_line(colour="black") +
  labs(x="Theoretical Cumulative Distribution",
       y="Empirical Cumulative Distribution") +
  theme(axis.title.x = element_text(size=6),
        axis.title.y = element_text(size=6),
        axis.text = element_text(size=5),
        plot.background = element_blank(),
        panel.background = element_rect(fill="white", colour="black"),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        text = element_text(family="Serif", face="bold"))

# append plots (create Figure S7)
fig_s7 <- ggarrange(pp_plot_m5_country1, pp_plot_m5_country2,
                    pp_plot_m5_country_year1, pp_plot_m5_country_year2,
                    ncol=2, nrow=2, labels=c("a","","b",""),
                    font.label=list(size=12, color="black"))

# save plot as .png file
ggsave(filename="Figures/fig_s7.png", plot=fig_s7, width=6, height=4, units="in", dpi=1500)







